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Summary. We investigate the stress response function of a layer of grains, i.e. 
the stress profile in response to a localized overload. The shape of the profile is very 
sensitive to the packing arrangement, and is thus a good signature of the preparation 
procedure of the layer. This study has been done by the use of molecular dynamics 
numerical simulations. Here, for a given rain-like preparation, we present the scaling 
properties of the response function, and in particular the influence of the thickness 
of the layer, and the importance of the location of the overload and measurement 
points (at the boundaries, in the bulk). 
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The statics of granular materials has been a rich field of research over the 
last few years. One of the reasons of this interest is that a system of grains 
reaches its mechanical equilibrium after the particles are 'jammed' in config- 
urations in relation to the previous dynamics of the grains. As a consequence, 
the distribution of the stresses in a static piling of grains depends, in a subtle 
manner, on the preparation procedure of the system. A now famous example 
is that of the sand pile for which the pressure profile at the bottom is differ- 
ent when built from a hopper, or by successive horizontal layers (a 'rain' of 
grains). In the former situation, this profile exhibits a clear 'dip' below the 
apex of the pile [2], whereas it shows a flat 'hump' in the latter [1]. 

In fact, more interesting than the case of a pile, is the study of the stress 
response function of a layer of grains, i.e. the stress profiles in response to a 
localized overload, see figure 1. This geometry allows much more variety in the 
preparation of the layer. Experiments have been performed for instance with 
compacted, loose or sheared layers [3-5], whose response pressure profiles are 
different enough from each other to be a kind of 'signature' of their texture. 
Besides, such a response test is the elementary 'brick' with which other sit- 
uations can be deduced, e.g. that of the pile [6]. It is then also well adapted 
for the comparison of the different models of stress distribution in granular 
media. 
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Fig. 1. System geometry and notations. The stress measure is done at a position 
(a;, z) from the applied force Fo and at a distance h — z from the bottom of the layer. 
To compute the stress, contact forces are integrated over a linear extension I. Note 
that we use horizontal periodic boundary conditions. 



Our aim is to perform extensive simulations of assemblies of grains, in order 
to provide precise two-dimensional numerical data of stress response functions. 
The control of all the parameters of the simulations, as well as the ability of 
measuring both micro (grain size) and macro (system size) quantities, ensure 
a useful and interesting feed back to the experiments and the models. 

Although our ambition is to be as general and systematic as possible, we 
shall in this paper, present only few results concerning some scaling and size 
effects of this response function. In particular we shall restrict to a single (rain- 
like) preparation history and study the influence of the thickness of the layer, 
as well as that of the location of the overload and measure points. Besides, 
the overloading procedure requires careful and important validity tests (e.g. 
linearity, additivity, reversibility) that will be also described below. 

1 The simulations 

In figure 2, we show a part of a typical simulated layer, where we can see the 
initial force chain network and the response produced by a localized overload 
applied to a single grain after the initial forces have been subtracted. The sim- 
ulations are performed using a classical molecular dynamics (MD) algorithm 
in three successive stages: preparation, deposition and overloading that are 
described in the following. 

We start with the preparation of N polydisperse grains, with radii homoge- 
neously distributed between R m in and R m ax = 2i? m i n . These N particles are 
put on the nodes of a grid with aspect ratio 1:4- this aspect ratio is needed 
for the appropriate study of the response functions, see below. This lattice 
ends at its bottom with a fixed horizontal line of similar particles which will 
be used as the support for deposition - the distances between these bottom 
grains are small enough to avoid grain evasion. 

The deposition stage consists of untying the particles from the grid and, un- 
der gravity and horizontal periodic boundary conditions, letting them evolve 
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Fig. 2. Example of the granular layers produced by our simulations. Left: Force 
chain network after the deposition under gravity of N = 1600 grains. This picture 
corresponds to the equilibrium configuration after 628100 MD steps. Right: Response 
due to an overload localized over a single grain after subtraction of the initial forces. 
The amplitude of this overload corresponds to the weight of five average grains. 
This picture has been obtained after 29500 additional MD steps to reach a new 
equilibrium configuration. In black, contact forces have been increased in response 
to the overload, whereas in gray they have been decreased. 

with random initial velocities. In this 'rain-like' deposition the grains interact 
with elastic and friction forces, and the system evolves under classical New- 
ton equations and molecular dynamics rules (e.g. predictor-corrector, Verlet 
algorithm). The time to reach the equilibrium depends on the characteristics 
of the particles. We chose /i = 0.5 for the contact friction coefficient. The 
rhcology of the contacts is that of Kelvin- Voigt with k n = 1000 N/m and 
k t = 750 N/m for the normal and horizontal contact stiffness. The viscosity 
g n is chosen in order to get a critical damping. At last, the gravity is set to 
unity. The equilibrium criteria consist of the five following tests which are 
applied after each period of 100 MD time steps: (1) the number of gained/lost 
contacts during this period has to be zero; (2) the number of sliding contacts 
between particles also has to vanish; (3) the integrated force measured at the 
bottom of the layer must be equal to the sum of the weight of all the grains; 
(4) all the particles have to have, at least, two contacts; and (5) the total ki- 
netic energy has to be lower than some low threshold. Once these criteria are 
all satisfied, the deposition is stopped, and the overloading phase can begin. 

The overload is applied to a single grain with a force F that we wish 
to be small enough to not cause any rearrangement of the layer structure. 
The determination of the amplitude of the loading force is then a crucial part 
of the simulation. For that purpose we made several tests which check the 
reversibility, the additivity and the linearity of the stress response. These tests 
(see below for more details) let to conclude that the optimal magnitude of Fq 
is of the order of few times the weight of an average grain. In the overloading 
procedure, Fq is split in 10000 MD time steps, i.e. in each time step the force 
is increased in 1/10000 of Fq, and it remains constant and equal to Fo after 
that. Again, the equilibrium criteria are those explicited above, except that 
the integrated force at the bottom of the layer has to be equal to the weight 
of all the grains plus F . 
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Fig. 3. In this picture we show one of the tests applied to the layers (linearity, 
additivity and reversibility) to settle the optimal value for Fo. We present here the 
linearity test which consists in comparing the force response profiles with different 
values of Fo. These curves have been obtained on a layer with N = 400 grains 
(h ~ ID), and F m — 2.5 average grain weight. Rescaled by the overload force, the 
difference of the two responses is reasonnably small. This difference is due to small 
slips in the contacts between the grains. 

Only the linearity test is illustrated here, in figure 3: it consists in compar- 
ing the force response after loading with different values of Fo. We find that the 
response is satisfactorily linear in F as long as its magnitude does not exceed 
~ 40 times the average grain weight. With the other tests we have checked 
that loading and unloading with the very same F gives the original contact 
force distribution with a good precision, and that two simultaneous overloads 
located at two different positions give a response which is perfectly equivalent 
to the sum of the two responses computed from the two corresponding single 
overloads. 



2 Results and discussion 

In order to measure the vertical stress a zz , we integrate the contact forces 
over a set of grains. In fact, we do not make use of the usual formula cr Q( g — 
l/SJ2i fa r pi where the sum is computed over the contacts between the grains 
in the 'volume' S (here we are two-dimensional), f a being a th component of 
the considered contact force, and rp the [3 th component of the corresponding 
distance vector between the grains in contact. Rather, we take a horizontal 




Fig. 4. Response functions. Several curves are shown for different values of z and h. 
On the left, the response profiles collapse when measured at the bottom of the layers 
and rescaled by h. On the right, the response in the bulk. Although more averaging 
is needed, this graph suggests interesting features for z < h. All these curves have 
been computed with i — 0.772. 



'line' of grains of length £. The corresponding stress a zz is then equal to the 
sum of the vertical components of the forces carried by the contacts of one of 
the sides of the line - say the upper one - divided by I. This stress measure can 
be done at any depth z (z = h means that the measure is done on the bottom) 
and centered at any horizontal distance x from the overload point. As we are 
particulary interested in a zz profiles along a: at a given z, our choice is here 
better adapted than the usual stress formula which would mix together grains 
of (slightly) different depths in S. This is particulary important at small z. 
Besides, with this stress definition, the integral of <r zz over x at a fixed z is 
exactly equal to Fq at any scale. This property is crucial to normalize and 
compare data from layers of different thickness. 

The response function at depth z is obtained by making the difference of 
the stress profile a zz (x, z) measured on the layer with the additional force Jo, 
with a copy of the very same layer without overload. This difference is ensem- 
ble averaged over several overload positions and layer samples. We present 
here the results for five different layer sizes, with N = 100,400,1600,3600 
and 6400 grains. The 1 : 4 aspect ratio gives layers of average thickness of 
4, 7, 14, 21 and 29 mean particle diameters D, respectively. The 10, 5, 5, 4 and 
4 different realizations give in the end 62, 60, 65, 84 and 86 different overload 
points - on average, there are six grains between two consecutive overload 
points. 

The results are shown on figure 4. One can see that the normalized pro- 
files measured at the bottom of layers of different thickness h collapse when 
they arc linearly rescaled by h, and show a single peak whose width is of the 
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order of h. This is precisely the reason why it is important to have layers four 
times wider than thick - 2h on both sides of the overload point. When the 
response is computed in the bulk, the structure of the response profile looks 
more complicated. The present curve needs however more ensemble averaging 
and should be taken as preliminary. These results were all obtained using a 
measure scale I = 0.77z for all layer sizes. A systematic study of the impor- 
tance of this integration measure length I is under way. We are also working 
on the calculation of the other stress components a xx and a xz ■ 

These scaling studies are currently extended to other preparation histories, 
and in particular to the cases of sheared or more anisotropic layers of grains. At 
last, we plan to test these data against the predictions of anisotropic elasticity 
which, depending on the values of the different parameters, can give various 
stress response profiles [7] . 
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